# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Collective victimhood beliefs and conflict-related attitudes: A meta-analysis
# Marko Kljajic, Nadav Shelef, and Ethan vanderWilden

# Last Replicated: August 6, 2025

# Analysis file
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

############ Load in packages, data, helper functions ############
# List of packages
pkg = c("tidyverse", "patchwork", "metafor", "kableExtra", "rmeta", "reshape2", "metaSEM")

# Checks if they are installed, install if not
if (length(setdiff(pkg, rownames(installed.packages()))) > 0) {install.packages(setdiff(pkg, rownames(installed.packages()))) }

# Load + clear environment
suppressMessages(suppressWarnings(lapply(pkg, library, character.only = TRUE)))
rm(list=ls())

# set working directory (NOTE: change for local file path)
setwd("C:/Users/19785/Desktop/CPS_replication")

# load in the dataset
data <- read.csv("KSV_dataset.csv")

#Load in functions from file
source("helperFunctions.R")

### Main analysis excludes 'bundled' treatments ; emotions/cognitive perspective outcomes
main <- subset(data, data$include.main == 1) # (keep emotions and cognitive frames)
main <- subset(main, main$outcome %in% c("Hawkishness", "Ingroup Attachment", "Outgroup Exclusion", "Reconciliation"))



############ Appendix Analyses ####

### FIGURE A1 - attention to CV , RENDER 12 x 6
data %>% group_by(Year) %>% dplyr::summarise(estimates = n()) -> time
data %>% group_by(Title, Year) %>% 
  dplyr::summarise(estimates = n()) %>% group_by(Year) %>%
  dplyr::summarise(manuscripts = n()) %>% merge(time, by = c("Year")) -> time

cols = c("Estimates" = "black", "Manuscripts" = "darkgrey")

ggplot(data = time) +
  geom_line(aes(x = Year, y = estimates, color = "Estimates"), linewidth = 1.2) +
  geom_line(aes(x = Year, y = manuscripts/0.1, color = "Manuscripts"), linewidth = 1.2) +
  scale_x_continuous(limits = c(2007, 2023), breaks = seq(2007, 2023, 2)) +
  scale_y_continuous(name = "Estimates", limits = c(0, 330), breaks = seq(0, 300, 100),
                     sec.axis = sec_axis(~.*0.1, name = "Manuscripts")) +
  scale_color_manual(values = cols) + theme_bw() +
  theme(legend.title = element_blank(),
        axis.text.x = element_text(angle = 45, vjust = 0.5, size = 12),
        axis.text.y = element_text(12),
        axis.title = element_text(size = 14),
        legend.position = c(0.2, 0.85),
        legend.text = element_text(size = 14)) 





### FIGURE A2 - Relaxed inclusion criteria + Remove Victim-Perpetrator comparisons , RENDER 20 x 10

#remove emotions and cognitive perspectives, but broaden inclusion criteria
robust <- subset(data, data$outcome %in% c("Hawkishness", "Ingroup Attachment", "Outgroup Exclusion", "Reconciliation"))

full_plots(fill_metatable3(robust, exclude_mechs = T), df = robust,
           exclude_mechs = T, Bonf = T, range = 1.5) +
  theme(legend.position = "none", plot.title = element_text(size = 16, face = 'bold', hjust = 0.5)) +
  labs(title = "Relaxed inclusion criteria") +
  
  #remove victim-perpetrator plots
  full_plots(fill_metatable3(subset(main, main$VP.comparison == 0), exclude_mechs = T),
             df = subset(main, main$VP.comparison == 0), exclude_mechs = T, Bonf = T, range = 1.5) +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), 
        plot.title = element_text(size = 16, face = 'bold', hjust = 0.5)) +
  labs(title = "Remove victim-perpetrator comparisons")







### FIGURE A3 Removing Outliers  , RENDER 20 x 10

# first get meta results
results <- fill_metatable3(main, exclude_mechs = T)

#Log whether an estimate overlaps with the corresponding overall victimhood CI
main$overlap <- rep(NA, nrow(main))

for (i in 1:nrow(main)){
  if (main$treatment[i] == "Collective victimhood"){
    est <- as.numeric(strsplit(results[2,1], " ")[[1]][1])
    se <- as.numeric(substr(strsplit(results[2,1], " ")[[1]][2], 2, 6))
    main$overlap[i] <- ifelse(main$d.directional[i] + 1.645*main$sd.d[i] < est - 1.645*se | 
                                main$d.directional[i] - 1.645*main$sd.d[i] > est + 1.645*se, 0, 1)}
  else if (main$treatment[i] == "Competitive victimhood"){
    est <- as.numeric(strsplit(results[3,1], " ")[[1]][1])
    se <- as.numeric(substr(strsplit(results[3,1], " ")[[1]][2], 2, 6))
    main$overlap[i] <- ifelse(main$d.directional[i] + 1.645*main$sd.d[i] < est - 1.645*se | 
                                main$d.directional[i] - 1.645*main$sd.d[i] > est + 1.645*se, 0, 1)}
  else if (main$treatment[i] == "Exclusive victimhood"){
    est <- as.numeric(strsplit(results[4,1], " ")[[1]][1])
    se <- as.numeric(substr(strsplit(results[4,1], " ")[[1]][2], 2, 6))
    main$overlap[i] <- ifelse(main$d.directional[i] + 1.645*main$sd.d[i] < est - 1.645*se | 
                                main$d.directional[i] - 1.645*main$sd.d[i] > est + 1.645*se, 0, 1)}
  else if (main$treatment[i] == "Inclusive victimhood"){
    est <- as.numeric(strsplit(results[5,1], " ")[[1]][1])
    se <- as.numeric(substr(strsplit(results[5,1], " ")[[1]][2], 2, 6))
    main$overlap[i] <- ifelse(main$d.directional[i] + 1.645*main$sd.d[i] < est - 1.645*se | 
                                main$d.directional[i] - 1.645*main$sd.d[i] > est + 1.645*se, 0, 1)}
}

#plot when removing 5% most extreme estimates
full_plots(fill_metatable3(subset(main, main$d < 1.3), exclude_mechs = T), #remove top 5 % (37 estimates)
           df = subset(main, main$d < 1.3), exclude_mechs = T, Bonf = T, range = 1.5) +
  theme(legend.position = "none", plot.title = element_text(size = 16, face = 'bold', hjust = 0.5)) +
  labs(title = "Remove most extreme 5% of estimates") +
  
  #plot when removing non-overlapping estimates
  full_plots(fill_metatable3(subset(main, main$overlap == 1), exclude_mechs = T), #remove non-overlapping estimates
             df = subset(main, main$overlap == 1), exclude_mechs = T, Bonf = T, range = 1.5) +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), 
        plot.title = element_text(size = 16, face = 'bold', hjust = 0.5)) +
  labs(title = "Remove non-overlapping estimates")






### FIGURE A4 - remove the US and Israel , RENDER 20 x 10
#Removes the US
full_plots(fill_metatable3(subset(main, main$Setting != "United States"), exclude_mechs = T), 
           df = subset(main, main$Setting != "United States"), exclude_mechs = T, Bonf = T, range = 1.5) +
  theme(legend.position = "none", plot.title = element_text(size = 16, face = 'bold', hjust = 0.5)) +
  labs(title = "Remove estimates from the United States") +
  
  #Remove Israel
  full_plots(fill_metatable3(subset(main, main$Setting != "Israel"), exclude_mechs = T),
             df = subset(main, main$Setting != "Israel"), exclude_mechs = T, Bonf = T, range = 1.5) +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), 
        plot.title = element_text(size = 16, face = 'bold', hjust = 0.5)) +
  labs(title = "Remove estimates from Israel")







### FIGURE A5 - Alternative modeling choices (NOTE: using functions 'fill_metatable' and 'fill_metatable2')
# Inverse variance weighted two-level model
full_plots(fillMeta_simple(main, exclude_mechs = T, weights = T, ivw = T, cap = 100), df = main, exclude_mechs = T, Bonf = T) +
  theme(legend.position = "none", plot.title = element_text(size = 16, face = 'bold', hjust = 0.5)) +
  labs(title = "Two-level model with inverse variance weighting") +
  
  # 2-stage meta model (Paluck et al. 2021)
  full_plots(fillMeta_twostage(main, exclude_mechs = T), df = main, exclude_mechs = T, Bonf = T) +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), 
        plot.title = element_text(size = 16, face = 'bold', hjust = 0.5)) +
  labs(title = "Two-stage meta analysis model")








### FIGURE A6 - Articles with pre-registrations , RENDER: 20 x 8
Edata <- subset(main, main$Design.binary == "Experimental")

nonprereg <- subset(Edata, Edata$PreRegistered == 0)
NPRoutputs <- fill_metatable3(nonprereg, exclude_mechs = T)

prereg <- subset(Edata, Edata$PreRegistered == 1)
PRoutputs <- fill_metatable3(prereg, exclude_mechs = T)

newer <- subset(Edata, Edata$Year >= 2020)
neweroutputs <- fill_metatable3(newer, exclude_mechs = T)

older <- subset(Edata, Edata$Year < 2020)
olderoutputs <- fill_metatable3(older, exclude_mechs = T)

comparison_plot(PRoutputs, NPRoutputs, prereg, nonprereg,
             "Pre-registered", "No Pre-registration",
             exclude_mechs = T, range = 0.75, jump = 0.25) + theme(legend.position = "top") +
  
  comparison_plot(olderoutputs, neweroutputs, older, newer,
               "Before 2020", "2020 or later",
               exclude_mechs = T, range = 0.75, jump = 0.25) +
  theme(legend.position = "top", axis.text.y = element_blank()) # 20 x 8






### FIGURE A7 - Regime type/democracy , RENDER 12 x 6
dems <- fill_metatable3(subset(main, main$democracy == 1), exclude_mechs = T)
nondems <- fill_metatable3(subset(main, main$democracy == 0), exclude_mechs = T)
comparison_plot(dems, nondems, subset(main, main$democracy == 1), subset(main, main$democracy == 0),
             "Democracy", "Non-Democracy",
             exclude_mechs = T, range = 1, jump = 0.25) + theme(legend.position = "top")






### FIGURE A7 - Original PAP figure , RENDER 14 x 14
full_plots(fillMeta_simple(data), df = data) 


